Contents

rm(list=ls())

## load libraries
library(readxl)
library(kableExtra)
library(ComplexHeatmap)
library(OlinkAnalyze)
library(factoextra)
library(ggpubr)
library(tidyverse)
library(RColorBrewer)
library(limma)
library(ggrepel)

filter = dplyr::filter
select = dplyr::select

palette_treatment <- c(brewer.pal(3, "Set1"))

col_treatment <- list(treatment = c("Healthy" = "white", "Vehicle" = "black", "Isotype" = "#FFD500", "aCD47" = "#F99C1C", "EGFRvIII CAR" = "#008FD5", "EGFRvIII CAR + aCD47" = "#009B49", "EGFRvIII-SGRP CAR" = "#0054A6", "CD19 CAR"= "#EE1C24", "CD19-SGRP CAR"= "#981B1E"))

col_treatment_dataset <- list(treatment = c("Healthy" = "white", "Vehicle" = "black", "Isotype" = "#FFD500", "aCD47" = "#F99C1C", "EGFRvIII CAR" = "#008FD5", "EGFRvIII CAR + aCD47" = "#009B49", "EGFRvIII-SGRP CAR" = "#0054A6", "CD19 CAR"= "#EE1C24", "CD19-SGRP CAR"= "#981B1E"), dataset = c("nov2022" = "lightgrey", "june2023" = "darkgrey", "bridged" = "#3f3f3f")) 

my_boxplot_colors <- scale_fill_manual(breaks = c("Healthy", "Vehicle", "Isotype", "aCD47", "EGFRvIII CAR", "EGFRvIII CAR + aCD47", "EGFRvIII-SGRP CAR", "CD19 CAR", "CD19-SGRP CAR"), 
                       values=c("white", "black", "#FFD500", "#F99C1C", "#008FD5", "#009B49", "#0054A6", "#EE1C24", "#981B1E"))

my_pca_colors <- scale_color_manual(breaks = c("Healthy", "Vehicle", "Isotype", "aCD47", "EGFRvIII CAR", "EGFRvIII CAR + aCD47", "EGFRvIII-SGRP CAR", "CD19 CAR", "CD19-SGRP CAR"), 
                       values=c("#dbdad7", "black", "#FFD500", "#F99C1C", "#008FD5", "#009B49", "#0054A6", "#EE1C24", "#981B1E"))

treatment_levels <- c("Healthy", "Vehicle", "Isotype", "aCD47", "CD19 CAR", "CD19-SGRP CAR", "EGFRvIII CAR", "EGFRvIII CAR + aCD47", "EGFRvIII-SGRP CAR")

## Setting up directories

dataDir <- "data"
rdsDir <- "rds"
plotDir <- "plots"

1 Introduction

human tumor cells injected in mice lacking Tcells, Bcells, NKcells, but have microglia and macrophages.

Plasma samples from immuno-compromised mice:

2 Dataset

my_NPX_data0A <- read_NPX(file.path(dataDir, "Nov2022 Immuno Oncology_Martins_NPX.xlsx"))
my_NPX_data0B <- read_NPX(file.path(dataDir, "Immuno Oncolog_Martins_NPX.xlsx"))


metadata0 <- read_xlsx(file.path(dataDir, "Sample list TM.xlsx"), sheet = 1)%>%
  dplyr::rename(SampleID=Name, treatment=Condition)%>%
  arrange(SampleID)%>%
  mutate(SampleID=gsub(" ","_",paste0("mouse_", SampleID), fixed=TRUE),
         treatment = factor(treatment, levels = treatment_levels))

my_NPX_dataA <- my_NPX_data0A%>%
  mutate(SampleID=gsub("B_","", gsub(" ","_",paste0("mouse_", SampleID), fixed=TRUE)))%>%
  left_join(metadata0, by="SampleID")%>%
  mutate(under_LOD=ifelse(NPX<LOD, "yes", "no"),
         dataset="nov2022",
         treatment = factor(treatment, levels = treatment_levels))

my_NPX_dataB <- my_NPX_data0B%>%
  mutate(SampleID=gsub("B_","", gsub(" ","_",paste0("mouse_", SampleID), fixed=TRUE)))%>%
  left_join(metadata0, by="SampleID")%>%
  mutate(under_LOD=ifelse(NPX<LOD, "yes", "no"),
         dataset="june2023",
         treatment = factor(treatment, levels = treatment_levels))

my_NPX_data <- bind_rows(my_NPX_dataA, my_NPX_dataB)

table(metadata0$SampleID %in% my_NPX_data$SampleID)
## 
## TRUE 
##   54

3 Plot NPX density before bridging normalization

my_NPX_data %>%
  mutate(Panel = gsub("Olink ", "", Panel)) %>%
  ggplot(aes(x = NPX, fill = dataset)) +
  geom_density(alpha = 0.4) +
  olink_fill_discrete(coloroption = c("red", "darkblue")) +
  set_plot_theme() +
  ggtitle("Before bridging normalization: NPX distribution") +
  theme(axis.title.x = element_blank(),
        axis.title.y = element_blank(),
        strip.text = element_text(size = 16),
        legend.title = element_blank(),
        legend.position = "top")

3.1 before bridging

#### Extract bridging samples

overlapping_samples <-  data.frame(SampleID = intersect(my_NPX_dataA$SampleID, my_NPX_dataB$SampleID)) %>%  
  filter(!str_detect(SampleID, "CONTROL_SAMPLE")) %>% #Remove control samples
  pull(SampleID)


npx_before_br <- my_NPX_dataA %>%
  dplyr::filter(!str_detect(SampleID, "CONTROL_SAMPLE")) %>% #Remove control samples
  dplyr::mutate(Type = if_else(SampleID %in% overlapping_samples,
                        paste0("nov2022 Bridge"),
                        paste0("nov2022 Sample"))) %>%
  rbind({
    my_NPX_dataB %>%
      filter(!str_detect(SampleID, "CONTROL_SAMPLE")) %>% #Remove control samples %>% 
      mutate(Type = if_else(SampleID %in% overlapping_samples,
                            paste0("june2023 Bridge"),
                            paste0("june2023 Sample"))) %>%
      mutate(SampleID = if_else(SampleID %in% overlapping_samples,
                                paste0(SampleID, "_new"),
                                SampleID))
  })

md_before_br <- npx_before_br%>%
  select(SampleID, treatment, dataset, Type, `Color code`)%>%
  unique()%>%
  arrange(SampleID)

3.2 PCA before bridging

NPX_noNA <-npx_before_br%>%
  dplyr::select(SampleID, Assay, NPX)%>%
  spread(Assay,NPX)%>%
  select_if(~ !any(is.na(.)))%>%
  arrange(SampleID)%>%
  column_to_rownames("SampleID") 

all.equal(rownames(NPX_noNA), md_before_br$SampleID)
## [1] TRUE
pca <- prcomp(NPX_noNA, scale = T)  

fviz_eig(pca) 

fviz_pca_ind(pca, axes = c(1, 2), label="none", habillage = md_before_br$Type,
                       addEllipses=TRUE, ellipse.level=0.95, pointsize = 2,
             invisible="quali")

3.2.1 PCA november 2022

NPX_noNA1 <-NPX_noNA[rownames(NPX_noNA) %in% filter(npx_before_br, dataset=="nov2022")$SampleID,]

md_before_br1 <- md_before_br%>%
  filter(dataset=="nov2022")

all.equal(rownames(NPX_noNA1), md_before_br1$SampleID)
## [1] TRUE
pca1 <- prcomp(NPX_noNA1, scale = T)  # x is a numeric matrix (row names = patient, colnames= molecule) 

fviz_eig(pca1) 

fviz_pca_ind(pca1, axes = c(1, 2), label="none", habillage = md_before_br1$Type,
                       addEllipses=TRUE, ellipse.level=0.95, pointsize = 2,
             invisible="quali")

3.2.2 PCA june 2023

NPX_noNA2 <-NPX_noNA[rownames(NPX_noNA) %in% filter(npx_before_br, dataset=="june2023")$SampleID,]

md_before_br2 <- md_before_br%>%
  filter(dataset=="june2023")

all.equal(rownames(NPX_noNA2), md_before_br2$SampleID)
## [1] TRUE
pca2 <- prcomp(NPX_noNA2, scale = T)  # x is a numeric matrix (row names = patient, colnames= molecule) 

fviz_eig(pca2) 

fviz_pca_ind(pca2, axes = c(1, 2), label="none", habillage = md_before_br2$Type,
                       addEllipses=TRUE, ellipse.level=0.95, pointsize = 2,
             invisible="quali")

3.3 distribution

npx_before_br %>%
  mutate(Panel = gsub("Olink ", "", Panel)) %>%
  ggplot(aes(x = NPX, group = SampleID, fill=dataset)) +
  geom_density(alpha = 0.4) +
  set_plot_theme() +
  ggtitle("Before bridging normalization: NPX distribution") +
  theme(axis.title.x = element_blank(),
        axis.title.y = element_blank(),
        strip.text = element_text(size = 16),
        legend.title = element_blank(),
        legend.position = "top")

4 Perform bridging between two data sets

reference dataset = november2022

overlap_samples_list <- list("DF1" = overlapping_samples,
                             "DF2" = overlapping_samples)

# Perform Bridging normalization
npx_br_data <- olink_normalization_bridge(project_1_df = my_NPX_dataA,
                                          project_2_df = my_NPX_dataB,
                                          bridge_samples = overlap_samples_list,
                                          project_1_name = "nov2022",
                                          project_2_name = "june2023",
                                          project_ref_name = "nov2022")
# dplyr::glimpse(npx_br_data)

npx_br_data %>%
  mutate(Panel = gsub("Olink ", "", Panel)) %>%
  ggplot2::ggplot(ggplot2::aes(x = NPX, fill = dataset)) +
  ggplot2::geom_density(alpha = 0.4) +
  olink_fill_discrete(coloroption = c("red", "darkblue")) +
  set_plot_theme() +
  ggplot2::ggtitle("After bridging normalization: NPX distribution") +
  theme(axis.title.x = element_blank(),
        axis.title.y = element_blank(),
        strip.text = element_text(size = 16),
        legend.title = element_blank(),
        legend.position = "top")

npx_after_br <- npx_br_data %>%
  dplyr::mutate(Type = ifelse(SampleID %in% overlapping_samples, 
                              paste(Project, "Bridge"),
                              paste(Project, "Sample")))%>%
  dplyr:::mutate(SampleID.orig = SampleID, dataset,
                 SampleID = paste0(SampleID, "_", dataset))


NPX_noNA <- npx_after_br%>%
  dplyr::select(SampleID, Assay, NPX)%>%
  spread(Assay, NPX)%>%
  select_if(~ !any(is.na(.)))%>%
  arrange(SampleID)%>%
  column_to_rownames("SampleID") 

md_after_br <- npx_after_br%>%
  select(SampleID, SampleID.orig, treatment, dataset, Type, `Color code`)%>%
  unique()%>%
  arrange(SampleID)

5 Whole dataset overview

5.1 Plot densities

NPX <-npx_after_br%>%
  dplyr::select(SampleID, Assay, NPX)%>%
  spread(Assay,NPX)%>%
  arrange(SampleID)%>%
  column_to_rownames(var = "SampleID") 

all.equal(rownames(NPX), md_after_br$SampleID)
## [1] TRUE
plotDensities(t(NPX), legend=F)

plotDensities(t(NPX), legend = "topright", group= md_after_br$treatment)

5.2 create unique data point per mouse

i.e: use mean of bridged samples

npx_after_br_mean <- npx_after_br%>%
  group_by(SampleID.orig, Assay)%>%
  mutate(NPX.ogig=NPX,
         NPX=mean(NPX))%>%
  ungroup()%>%
  distinct(SampleID.orig, Assay, .keep_all = TRUE)%>%
  mutate(SampleID=ifelse(SampleID.orig %in% overlapping_samples, SampleID.orig, SampleID),
         dataset=ifelse(SampleID.orig %in% overlapping_samples, "bridged", dataset))

NPX_mean <-npx_after_br_mean%>%
  dplyr::select(SampleID, Assay, NPX)%>%
  spread(Assay,NPX)%>%
  arrange(SampleID)%>%
  column_to_rownames(var = "SampleID") 

md_after_br_mean <- md_after_br%>%
  mutate(SampleID=ifelse(SampleID.orig %in% overlapping_samples, SampleID.orig, SampleID),
         dataset=ifelse(SampleID.orig %in% overlapping_samples, "bridged", dataset))%>%
  distinct(SampleID.orig, .keep_all = TRUE)

plotDensities(t(NPX_mean), legend=F)

plotDensities(t(NPX_mean), legend = "topright", group= md_after_br_mean$treatment)

6 Normalisation

6.1 using cyclic loess (as suggested by Julien)

NPX.norm_mean <- t(normalizeBetweenArrays(t(NPX_mean), method = "cyclicloess"))

all.equal(rownames(NPX.norm_mean), md_after_br_mean$SampleID)
## [1] TRUE
plotDensities(t(NPX.norm_mean), legend=F)

plotDensities(t(NPX.norm_mean), legend = "topright", group= md_after_br_mean$treatment)

NPX.norm_mean.lf <- data.frame(NPX.norm_mean, check.names = F)%>%
  rownames_to_column("SampleID")%>%
  pivot_longer(values_to = "NPX.norm", names_to = "Assay", !c(SampleID))

npx_after_br_mean <- npx_after_br_mean%>%
  left_join(NPX.norm_mean.lf, by=c("SampleID", "Assay"))%>%
  relocate(NPX.norm, .after=NPX)

rm(NPX.norm_mean.lf)

6.1.1 Heatmap showing all genes after normalisation

# annotation
lbls_col = gsub("_2$","", md_after_br_mean$SampleID.orig)
column_ha = HeatmapAnnotation(
  treatment = md_after_br_mean$treatment,
  dataset = md_after_br_mean$dataset,
  col = col_treatment_dataset, annotation_name_side="left")

## annotation missing freq
miss.per.dataset_mean<-npx_after_br_mean%>%
  group_by(dataset, Assay)%>%
  summarise(MissingFreq=MissingFreq)%>%
  unique()%>%
  pivot_wider(names_from = dataset, values_from = MissingFreq)%>%
  type.convert()


all.equal(miss.per.dataset_mean$Assay, colnames(NPX.norm_mean))
## [1] "5 string mismatches"
miss.per.dataset_mean<-miss.per.dataset_mean[match(colnames(NPX.norm_mean), miss.per.dataset_mean$Assay),]

all.equal(miss.per.dataset_mean$Assay, colnames(NPX.norm_mean))
## [1] TRUE
row_ha = rowAnnotation(
  MissingFreq_june2023 = anno_barplot(miss.per.dataset_mean$june2023, gp = gpar(fill = col_treatment_dataset$dataset["june2023"])), 
  MissingFreq_nov2022 = anno_barplot(miss.per.dataset_mean$nov2022, gp = gpar(fill = col_treatment_dataset$dataset["nov2022"])), annotation_name_side="top")

# heatmap
hm <- Heatmap(t(scale(NPX.norm_mean, center = T)),
               column_labels = lbls_col,
               bottom_annotation = column_ha,
               right_annotation = row_ha,
               na_col = "black",
               # col = col_fun,
               name = "Z-Score",
               column_names_gp = grid::gpar(fontsize = 10),
               row_names_gp = grid::gpar(fontsize = 7))
hm

pdf(file.path(plotDir, "Supp_Fig5d_Normalized_IO_prot_in_plasma.pdf"), height =10, width=10)
hm
dev.off()
## quartz_off_screen 
##                 2

6.1.2 PCA after normalisation

remove NAs

NPX.norm_noNA_mean <- as.data.frame(NPX.norm_mean)%>%
  select_if(~ !any(is.na(.)))

all.equal(rownames(NPX.norm_noNA_mean), md_after_br_mean$SampleID)
## [1] TRUE
pca <- prcomp(NPX.norm_noNA_mean, scale = T) 
summary(pca)  
## Importance of components:
##                           PC1    PC2     PC3     PC4     PC5    PC6     PC7
## Standard deviation     3.9311 3.1320 2.55332 2.10463 1.96399 1.8687 1.72941
## Proportion of Variance 0.1717 0.1090 0.07244 0.04922 0.04286 0.0388 0.03323
## Cumulative Proportion  0.1717 0.2807 0.35313 0.40235 0.44521 0.4840 0.51724
##                            PC8     PC9    PC10    PC11    PC12    PC13    PC14
## Standard deviation     1.69400 1.68393 1.62095 1.55213 1.52804 1.46194 1.43134
## Proportion of Variance 0.03188 0.03151 0.02919 0.02677 0.02594 0.02375 0.02276
## Cumulative Proportion  0.54912 0.58063 0.60982 0.63659 0.66254 0.68628 0.70905
##                           PC15    PC16    PC17    PC18    PC19   PC20    PC21
## Standard deviation     1.36434 1.32921 1.28198 1.26903 1.25160 1.1503 1.14732
## Proportion of Variance 0.02068 0.01963 0.01826 0.01789 0.01741 0.0147 0.01463
## Cumulative Proportion  0.72973 0.74936 0.76762 0.78551 0.80292 0.8176 0.83225
##                           PC22    PC23    PC24    PC25    PC26    PC27    PC28
## Standard deviation     1.10338 1.05048 1.02724 0.97778 0.93495 0.92772 0.89838
## Proportion of Variance 0.01353 0.01226 0.01172 0.01062 0.00971 0.00956 0.00897
## Cumulative Proportion  0.84578 0.85804 0.86976 0.88038 0.89010 0.89966 0.90863
##                           PC29    PC30    PC31    PC32    PC33    PC34    PC35
## Standard deviation     0.87035 0.85432 0.81265 0.80314 0.76918 0.71921 0.71486
## Proportion of Variance 0.00842 0.00811 0.00734 0.00717 0.00657 0.00575 0.00568
## Cumulative Proportion  0.91704 0.92515 0.93249 0.93966 0.94623 0.95198 0.95766
##                           PC36    PC37    PC38   PC39    PC40    PC41    PC42
## Standard deviation     0.66055 0.64295 0.60276 0.5771 0.54667 0.53153 0.52491
## Proportion of Variance 0.00485 0.00459 0.00404 0.0037 0.00332 0.00314 0.00306
## Cumulative Proportion  0.96251 0.96710 0.97114 0.9748 0.97816 0.98130 0.98436
##                           PC43    PC44    PC45    PC46    PC47    PC48    PC49
## Standard deviation     0.48292 0.45429 0.42070 0.40445 0.37370 0.33992 0.31297
## Proportion of Variance 0.00259 0.00229 0.00197 0.00182 0.00155 0.00128 0.00109
## Cumulative Proportion  0.98695 0.98924 0.99121 0.99303 0.99458 0.99586 0.99695
##                           PC50    PC51    PC52    PC53      PC54
## Standard deviation     0.30463 0.28030 0.23838 0.21528 7.382e-15
## Proportion of Variance 0.00103 0.00087 0.00063 0.00051 0.000e+00
## Cumulative Proportion  0.99798 0.99885 0.99949 1.00000 1.000e+00
fviz_eig(pca) 

fviz_pca_ind(pca, axes = c(1, 2), label="none", habillage = md_after_br_mean$treatment,
                       addEllipses=TRUE, ellipse.level=0.95, pointsize = 2, invisible="quali")+
  my_pca_colors+
  my_boxplot_colors

pca1 <- as.data.frame(pca$x)%>%
  rownames_to_column("SampleID")%>%
  left_join(md_after_br_mean)

pdf(file.path(plotDir, "Fig4a_PCA_plasma_markers_distrib.pdf"), width=10, height = 6)
fviz_pca_ind(pca, axes = c(1, 2), label="none", habillage = md_after_br_mean$treatment,
                       addEllipses=TRUE, ellipse.level=0.80, pointsize = 2, invisible="quali")+
  my_pca_colors+
  my_boxplot_colors

ggplot(pca1, aes(PC1,PC2))+
  theme_bw()+
  geom_point(aes(color=treatment, shape=treatment), size=3)+
  stat_ellipse(aes(color=treatment, fill=treatment), geom = "polygon", alpha=0.2)+
  geom_hline(yintercept=0, linetype="dashed")+
  geom_vline(xintercept=0, linetype="dashed")+
  my_pca_colors+
  my_boxplot_colors+
  scale_shape_manual(breaks = c("Healthy", "Vehicle", "Isotype", "aCD47", "EGFRvIII CAR", "EGFRvIII CAR + aCD47", "EGFRvIII-SGRP CAR", "CD19 CAR", "CD19-SGRP CAR"), values=c(16, 17, 18, 0, 6, 8, 5, 7, 15))

dev.off()
## quartz_off_screen 
##                 2

7 Differential expression analysis

all.equal(md_after_br_mean$SampleID, rownames(NPX.norm_mean))
## [1] TRUE
moma <- model.matrix(~ 0 + treatment + dataset, data=md_after_br_mean)

colnames(moma) <- gsub("[^[:alnum:]]+", ".", gsub("treatment", "", colnames(moma), fixed = T))

## Contrast matrix
contrasts.matrix <- makeContrasts(
  EGFRvIII.CAR_vs_EGFRvIII.CAR.aCD47 =     EGFRvIII.CAR - EGFRvIII.CAR.aCD47,
  EGFRvIII.CAR_vs_EGFRvIII.SGRP.CAR =      EGFRvIII.CAR - EGFRvIII.SGRP.CAR,
  # EGFRvIII.CAR_vs_CD19.SGRP.CAR =          EGFRvIII.CAR - CD19.SGRP.CAR,
  EGFRvIII.CAR_vs_CD19.CAR =               EGFRvIII.CAR - CD19.CAR,
  CD19.SGRP.CAR_vs_EGFRvIII.SGRP.CAR =     CD19.SGRP.CAR - EGFRvIII.SGRP.CAR,
  CD19.CAR_vs_CD19.SGRP.CAR =              CD19.CAR - CD19.SGRP.CAR,
  EGFRvIII.CAR.aCD47_vs_EGFRvIII.SGRP.CAR =      EGFRvIII.CAR.aCD47 - EGFRvIII.SGRP.CAR,
  aCD47_vs_EGFRvIII.CAR.aCD47 =      aCD47 - EGFRvIII.CAR.aCD47,
  levels=moma)


## limma-trend
fit <- lmFit(t(NPX.norm_mean), moma)
fit2 <- contrasts.fit(fit, contrasts.matrix)
fit2 <- eBayes(fit2, trend=T, robust=F) ## robust=TRUE avoids shrinkage of large variances
table(de <- decideTests(fit2, p.value = 0.1))
## 
##  -1   0   1 
##  41 572  31

7.0.1 EGFRvIII vs EGFRvIII-SGRP

results <- topTable(fit2, p.value = 1, number = Inf, coef="EGFRvIII.CAR_vs_EGFRvIII.SGRP.CAR")

kable(results) %>% 
  kable_styling(position = "center", bootstrap_options = c("striped", "hover")) %>% 
  scroll_box(width = "800px", height = "400px")
logFC AveExpr t P.Value adj.P.Val B
CCL3 -0.9623920 -0.9667431 -6.6300760 0.0000000 0.0000028 8.7266751
CD27 0.5989700 0.0434190 4.2629293 0.0000968 0.0044517 0.8723171
IL13 -1.2176303 -0.7713638 -3.5113456 0.0009978 0.0305987 -1.3492629
FGF2 -0.3815010 -1.4941712 -2.8468567 0.0065331 0.1242215 -3.0987850
IL-1 alpha -0.7644022 -2.8846819 -2.8345475 0.0067512 0.1242215 -3.1288765
HGF 0.5530473 -0.6376821 2.5662891 0.0135365 0.2002119 -3.7611730
CD5 0.5656122 -0.2037086 2.5191764 0.0152335 0.2002119 -3.8674044
MUC-16 0.3009491 -1.6923834 2.4529876 0.0179448 0.2063649 -4.0141162
ARG1 -0.7029152 0.2162739 -2.3758325 0.0216507 0.2213182 -4.1813242
KLRD1 0.3425087 -0.8553006 2.2409556 0.0298070 0.2504705 -4.4634884
IL8 0.8367005 1.1624727 2.2389350 0.0299476 0.2504705 -4.4676154
CD40 0.5446794 -1.1657998 2.1036154 0.0408073 0.3128557 -4.7370979
IL18 -0.3684830 -1.6520022 -2.0368869 0.0473329 0.3349713 -4.8648831
FASLG -0.2419810 -0.7954715 -1.9414224 0.0582340 0.3652762 -5.0416837
HO-1 -0.8897826 1.8787318 -1.8997029 0.0636364 0.3652762 -5.1166867
MCP-2 0.3254890 -0.6541967 1.8677018 0.0680648 0.3652762 -5.1732729
IL12RB1 0.3427914 -1.2190962 1.8417367 0.0718477 0.3652762 -5.2185780
TRAIL -0.2379177 -1.0874565 -1.8401513 0.0720843 0.3652762 -5.2213265
CD8A 0.3620962 -1.4420661 1.8181448 0.0754375 0.3652762 -5.2592668
PTN 0.4286965 -2.2701437 1.7609400 0.0847733 0.3899571 -5.3560339
CASP-8 -0.5317449 -1.3645874 -1.6785476 0.0998980 0.4373142 -5.4906354
TWEAK 0.3084861 1.9511467 1.6551612 0.1045751 0.4373142 -5.5278021
TNF -0.2953175 -0.7343338 -1.5749503 0.1219961 0.4453495 -5.6517413
IL15 0.6558730 -0.7330178 1.5524183 0.1277490 0.4453495 -5.6557620
MMP7 0.6150454 1.6256308 1.5454253 0.1289709 0.4453495 -5.6959712
KIR3DL1 0.3366558 -2.5137428 1.5405396 0.1301554 0.4453495 -5.7032174
CRTAM -0.2660694 -1.2714628 -1.5383033 0.1307004 0.4453495 -5.7065273
Gal-9 0.1878269 -1.1532767 1.4830529 0.1447554 0.4756248 -5.7869183
TNFRSF21 0.2100938 5.6156598 1.4312399 0.1589974 0.5044055 -5.8598749
IL10 0.2432414 -1.5527216 1.3751041 0.1756385 0.5336925 -5.9362337
MMP12 -0.1643061 -1.3175578 -1.3616129 0.1798312 0.5336925 -5.9541659
CD4 0.2562845 -1.9008701 1.3349838 0.1883315 0.5414531 -5.9890811
PDGF subunit B -0.6601414 6.5580033 -1.3105465 0.1963983 0.5452567 -6.0205603
LAP TGF-beta-1 -0.5032287 3.4765154 -1.2875889 0.2042122 0.5452567 -6.0496413
MCP-3 0.2168590 -1.7680756 1.2783134 0.2074346 0.5452567 -6.0612551
ICOSLG -0.1447478 -1.3615353 -1.2372297 0.2221661 0.5677579 -6.1117531
NOS3 -0.3621967 -1.5350210 -1.1761168 0.2454888 0.6104047 -6.1840076
CCL23 -0.1538551 -0.4200597 -1.1272682 0.2653685 0.6424710 -6.2392806
MIC-A/B 0.1944967 -1.3480782 1.0539194 0.2973240 0.6882977 -6.3180982
DCN 0.1432047 -2.1245518 1.0324737 0.3071498 0.6882977 -6.3401889
CXCL13 -0.2062826 -2.0376865 -1.0183412 0.3137450 0.6882977 -6.3545094
VEGFA 0.1707794 -1.4025570 1.0173249 0.3142229 0.6882977 -6.3555319
IL33 0.1264109 -2.2348357 0.9458131 0.3490957 0.7302897 -6.4250240
CXCL10 -0.2895126 -2.1843327 -0.9318619 0.3561845 0.7302897 -6.4380139
IFN-gamma -0.7195532 4.2931944 -0.9193190 0.3627392 0.7302897 -6.4419794
CXCL12 0.1666708 -0.3822878 0.9144832 0.3651448 0.7302897 -6.4539349
ANGPT2 -0.1222052 -1.9008921 -0.8919054 0.3770009 0.7379592 -6.4741869
GZMH -0.1335740 -1.6666856 -0.7576131 0.4524758 0.8530743 -6.5844914
CCL19 -0.1248333 -2.0365889 -0.7510865 0.4563561 0.8530743 -6.5894066
GZMA 0.1848707 -2.2433253 0.7389426 0.4636273 0.8530743 -6.5984417
GZMB -0.0884843 -1.4916136 -0.6530481 0.5169167 0.8790867 -6.6582309
CCL17 0.2417116 1.1842148 0.6527681 0.5170957 0.8790867 -6.6584140
TNFSF14 0.1735614 0.6228065 0.6264163 0.5340818 0.8790867 -6.6752981
IL7 0.1096035 -0.0869920 0.6234895 0.5359862 0.8790867 -6.6771311
PDCD1 0.1895465 -0.8762071 0.5924066 0.5564263 0.8790867 -6.6960772
TIE2 0.1008270 -1.8918147 0.5780468 0.5660002 0.8790867 -6.7045077
CCL4 -0.0872828 -0.5336423 -0.5757577 0.5675339 0.8790867 -6.7058328
IL12 0.1883311 -1.6968559 0.5579164 0.5795574 0.8790867 -6.7159828
IL5 -0.1798104 -1.4721251 -0.5523870 0.5833086 0.8790867 -6.7190646
CSF-1 -0.1548806 -1.3511787 -0.5491891 0.5854834 0.8790867 -6.7208330
CAIX -0.0781767 -1.1660219 -0.5477622 0.5864551 0.8790867 -6.7216189
CD40-L 0.1030995 -2.1076827 0.4995181 0.6197518 0.8790867 -6.7469989
TNFRSF9 -0.1107819 -0.2974723 -0.4969672 0.6215359 0.8790867 -6.7482765
Gal-1 -0.1074847 0.2941486 -0.4936574 0.6238541 0.8790867 -6.7499246
CXCL1 0.1070344 -0.6135119 0.4870278 0.6285091 0.8790867 -6.7531928
PGF 0.0868608 -0.6154360 0.4785671 0.6344721 0.8790867 -6.7573003
TNFRSF4 -0.1413895 -1.9483573 -0.4530296 0.6526182 0.8790867 -6.7692657
CD83 0.0599345 -0.7571030 0.4504176 0.6544865 0.8790867 -6.7704529
CXCL5 -0.0559320 -2.1444526 -0.4436813 0.6593150 0.8790867 -6.7734833
ANGPT1 0.0500481 8.2381087 0.4180397 0.6778277 0.8908592 -6.7846038
LAMP3 0.0666158 -1.2008492 0.4047024 0.6875378 0.8908940 -6.7901284
CD28 0.0575505 -1.0302463 0.3671227 0.7151803 0.9138414 -6.8047376
CCL20 0.0732618 -1.1730268 0.3229705 0.7481544 0.9270812 -6.8200942
IL6 0.1885328 -0.3982624 0.3207243 0.7498453 0.9270812 -6.8208231
EGF 0.1112740 -0.1910594 0.3128634 0.7557727 0.9270812 -6.8233344
ADA 0.0582176 -2.2622091 0.2921652 0.7714500 0.9338605 -6.8296498
IL4 0.0915919 -1.1894395 0.2591126 0.7966833 0.9395703 -6.8388419
LAG3 0.0509785 -2.3299600 0.2392899 0.8119243 0.9395703 -6.8438274
NCR1 0.0453799 -1.6939142 0.2274449 0.8210673 0.9395703 -6.8466177
CXCL11 0.1011961 -0.4275654 0.2181226 0.8282807 0.9395703 -6.8487143
MCP-1 0.0870971 0.7340598 0.2158379 0.8300509 0.9395703 -6.8492148
MCP-4 -0.0378219 0.9295177 -0.2022974 0.8405600 0.9395703 -6.8520730
IL2 -0.0679478 -1.8765798 -0.1917905 0.8487350 0.9395703 -6.8541636
CD70 -0.0217853 -0.3738565 -0.1800771 0.8578686 0.9395703 -6.8563629
PD-L1 -0.0669698 -1.3769988 -0.1578187 0.8752778 0.9449407 -6.8601610
VEGFR-2 -0.0216313 -1.4390177 -0.1423947 0.8873788 0.9449407 -6.8624996
CD244 -0.0191346 -1.9857875 -0.1290871 0.8978412 0.9449407 -6.8643244
PD-L2 0.0112252 -1.8406661 0.1153445 0.9086649 0.9449407 -6.8660213
CX3CL1 -0.0138398 -1.3156902 -0.0984455 0.9219985 0.9449407 -6.8678464
CXCL9 0.0107076 -1.3096721 0.0954068 0.9243986 0.9449407 -6.8681440
ADGRG1 0.0117546 -1.5616801 0.0696628 0.9447585 0.9551405 -6.8702910
TNFRSF12A 0.0098702 4.7079791 0.0386012 0.9693723 0.9693723 -6.8719901 
# Categorize Results based on P-value & FDR for plotting
results$Color <- "NS or FC < 0.5"
results$Color[results$P.Value < 0.05] <- "P < 0.05"
results$Color[results$adj.P.Val < 0.05] <- "adj.P.Val < 0.05"
results$Color[results$adj.P.Val < 0.001] <- "adj.P.Val < 0.001"
results$Color[abs(results$logFC) < 0.5] <- "NS or FC < 0.5"
results$Color <- factor(results$Color,
                        levels = c("NS or FC < 0.5", "P < 0.05",
                                   "adj.P.Val < 0.05", "adj.P.Val < 0.001"))

p1 <- ggplot(results,
       aes(x = logFC, y = -log10(P.Value),
           color = Color, label = rownames(results))) +
    geom_vline(xintercept = c(0.5, -0.5), lty = "dashed") +
    geom_hline(yintercept = -log10(0.05), lty = "dashed") +
    geom_point() +
    labs(x = "Enriched in EGFRvIII-SGRP <- log2(FC) -> Enriched in EGFRvIII",
         y = "Significance, -log10(P)",
         color = "Significance") +
    scale_color_manual(values = c(`adj.P.Val < 0.001` = "dodgerblue",
                                  `adj.P.Val < 0.05` = "lightblue",
                                  `P < 0.05` = "orange2",
                                  `NS or FC < 0.5` = "gray"),
                       guide = guide_legend(override.aes = list(size = 4))) +
    scale_y_continuous(expand = expansion(mult = c(0,0.05))) +
    geom_text_repel(min.segment.length = 0.3) +
    theme_bw(base_size = 16) +
    theme(legend.position = "bottom")


p1

pdf(file.path(plotDir, "Fig4b_volcano_EGFRvIII_CAR_vs_EGFRvIII_SGRP_CAR.pdf"))
p1
dev.off()
## quartz_off_screen 
##                 2

8 signifiantly DE prot

interest.prot2 <- unique(rownames(topTable(fit2, p.value = 0.05, number = Inf, coef=c(1:7))))

comparisons <- list(c("EGFRvIII CAR", "EGFRvIII CAR + aCD47"), 
                       c("EGFRvIII CAR", "EGFRvIII-SGRP CAR"), 
                       c("EGFRvIII CAR", "CD19 CAR"), 
                       c("CD19-SGRP CAR", "EGFRvIII-SGRP CAR"), 
                       c("CD19 CAR", "CD19-SGRP CAR"), 
                       c("EGFRvIII CAR + aCD47", "EGFRvIII-SGRP CAR"), 
                       c("aCD47", "EGFRvIII CAR + aCD47"))

data <- npx_after_br_mean%>%
  filter(Assay %in% interest.prot2,
         !is.na(NPX))%>%
  mutate(Assay=factor(Assay, levels = interest.prot2),
         treatment=factor(treatment, levels = treatment_levels))

stat.test <- data%>%
  group_by(Assay)%>%
  rstatix::wilcox_test(NPX.norm~treatment,
                       comparisons = comparisons)%>%
  rstatix::add_xy_position()%>%
  mutate(grp=paste0(group1, "_", group2))%>%
  group_by(grp)%>%
  mutate(p.adj=p.adjust(p, method = "BH", n = length(p)),
         p.adj.signif=ifelse(p.adj<0.001, "***", ifelse(p.adj<0.01, "**", ifelse(p.adj<=0.05, "*", "ns"))))%>%
  ungroup()
  

write.csv(select(stat.test, -groups), file.path(plotDir,"stat_test.csv"), row.names = F )

p2 <- ggplot(data)+
  theme_bw()+
  geom_boxplot(aes(x=treatment, y=NPX.norm, fill=treatment), outlier.shape = NA, alpha=0.7)+
  my_boxplot_colors+
  geom_jitter(aes(x=treatment, y=NPX.norm, color=under_LOD), alpha=0.8)+
  scale_color_manual(values=c("black", "#999999"))+
  theme(axis.text.x=element_blank(), #remove x axis labels
        axis.ticks.x=element_blank())+ #remove x axis ticks
  ggtitle("normalised values")+
  facet_wrap(~Assay, scales="free", ncol=3)

p2

pdf(file.path(plotDir, "Fig4c_norm_exp_DE_plasma_markers.pdf"), width = 11, height = 13)
p2
dev.off()
## quartz_off_screen 
##                 2

9 IL6 plasma expression at d15

p3 <- ggplot(filter(npx_after_br_mean, Assay=="IL6"))+
  theme_bw()+
  geom_boxplot(aes(x=treatment, y=NPX.norm, fill=treatment), outlier.shape = NA, alpha=0.7)+
  my_boxplot_colors+
  geom_jitter(aes(x=treatment, y=NPX.norm, color=under_LOD), alpha=0.8)+
  scale_color_manual(values=c("black", "#999999"))+
  theme(axis.text.x=element_blank(), #remove x axis labels
        axis.ticks.x=element_blank())+ #remove x axis ticks
  ggtitle("IL6 normalised values")

p3

pdf(file.path(plotDir, "Supp_Fig7c_Plasma_IL6.pdf"))
p3
dev.off()
## quartz_off_screen 
##                 2